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Abstract 



We report numerical results for velocity correlations in dense, gravity-driven granular flow down 



an inclined plane. For the grains on the surface layer, our results are consistent with experimental 
measurements reported by Pouliquen. We show that the correlation structure within planes parallel 
to the surface persists in the bulk. The two-point velocity correlation function exhibits exponential 
decay for small to intermediate values of the separation between spheres. The correlation lengths 
identified by exponential fits to the data show nontrivial dependence on the averaging time At 
used to determine grain velocities. We discuss the correlation length dependence on averaging time, 



incline angle, pile height, depth of the layer, system size and grain stiffness, and relate the results to 
other length scales associated with the rheology of the system. We find that correlation lengths are 
typically quite small, of the order of a particle diameter, and increase approximately logarithmically 
with a minimum pile height for which flow is possible, h stop , contrary to the theoretical expectation 

O ' 

O . of a proportional relationship between the two length scales. 
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I. INTRODUCTION 

Dense, gravity- driven granular flows down an inclined plane achieved the 

status of a key model system because of their relevance to many geological and industrial 
applications. Provided that the surface of the incline is sufficiently rough and the flow height 
is small compared to the width and length of the chute, such that transients and side-wall 
effects can be neglected, the flow behavior is controlled by the angle of inclination, 9, and 
granular layer thickness h. For a given 9 the flow is possible only for h > h stop (9). Here 
^stop is often referred to as the deposit function, because it is approximately equal to the 
thickness of the deposit remaining on the plane once the flow is stopped either due to the 
decrease of 9 or decrease of h. Recent experimental and numerical studies confirmed that for 
spherical grains, the depth-averaged steady state velocity of the flow u can be related to the 
deposit function through the following relationship, originally proposed by Pouliquen 

u /HA*, (i) 



\/gh h stop (9) 

where g is the gravitational acceleration and (3 ~ 0.13. This relationship suggests that a 
single scaling length controls both the deposit function and the rheology. 

A potential explanation of this scaling length as arising from correlations in granular mo- 
tion was advanced by some of us p , linking a correlation length in the flow to a rheologically 
defined "viscosity length scale", l v , defined by the Bagnold scaling form [6( 

o*z = pili\ (2) 

for a granular system with bulk density p, flowing under a shear stress a xz at a shear rate 7. 
Alternatively, it has recently been arg 

of diameter d and grain density p g flowing while subject to a pressure P, the rheology is 
controlled by the local scaling variable 

which is a generalization of Eq. ([T]) . For steady state chute flow, these two quantities are 
related through 

rl Ua-nff 

(4) 




whre is the volume fraction, thus they can be used interchangeably in this case. The fact 
that the same rheology can be advanced with or without an assumed fundamental length 
scale (a "correlation length") raises questions about the correct interpretation of this length 
scale, and the nature of its relationship to flow-induced structures, if any Pouliquen p] has 
recently reported experimental measurements of two point velocity correlation functions at 
the flow surface, and obtained results supporting a possible connection between the observed 
correlation length and the deposit function. 

In this work, we repeat and extend the correlation analysis presented by Pouliquen [l| 
to layers that are far from the bottom and surface of the flow, where true bulk rheology 
is observed, by taking advantage of information available from discrete element (DEM) 
simulations that is difficult to obtain experimentally. We are able to reproduce and further 
analyze the experimental results, and ultimately make the following observations: 

• All two-point velocity correlation functions exhibit exponential decay with relative 
distance, generally with different values of the correlation length \ a p for different 
components (3 of the velocity and a of the displacement vector. 

• The measured correlation lengths \ a p vary with the velocity measurement time At 
used to calculate the velocities from displacements. In order to be able to compare 
correlation lengths at different depths and between different runs, it is necessary to 
adjust the measurement time such that the particles that are being measured expe- 
rience a fixed, predetermined strain e = 7AI This yields correlation lengths that 
are independent of the position of the measured layer in the bulk, or the overall flow 
height. 

• Correlation lengths in the bulk are uniquely determined by the incline angle 6 (or 
equivalently by h st0 p or I). They increase with decreasing 9. 

• Correlation lengths are typically quite small, of the order of a particle diameter, and 
increase approximately logarithmically with h stop , contrary to the theoretical expecta- 
tion of a proportional relationship between the two length scales. 

• For very high piles and angles close to the angle of repose, the flows excite a low fre- 
quency "breathing mode" with coherent motion normal to the surface layer, previously 
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observed by Silbert [10]. This limits our ability to probe the limit of large h stop , small 
/ in this geometry. 

The paper is organized as follows: In Section |H] we outline the simulation method, which 
has been described in detail elsewhere 0, and explain the analysis used to obtain the 
reported results. In Section fill Al we present the results of velocity correlations analysis in 
the surface layer and compare these results with those of Pouliquen In Section UlIBI we 
discuss the differences between correlations on the surface and in the bulk. We then report 
the results for velocity correlations in the bulk in Section fill CI We also discuss the velocity 
fluctuations in Appendix [X] and the effect of system size and stiffness of particles on the 
correlation length in Appendix El Finally, in Section HVl we present our conclusions. 

II. SIMULATION AND ANALYSIS METHOD 




FIG. 1: (Color online) Typical snapshots from simulations with labeled coordinate axes. 9 is an 
angle of inclination of the rough base. Base area 20<i x 20d. 9089 spheres in the system of height 
h/d ~ 20 (left) and 17889 spheres in the system of height h/d ~ 40 (right). 



A. Simulation Method 

We perform three-dimensional discrete element (DEM) simulations of spherical, monodis- 
perse particles of diameter d in the chute flow geometry. The simulation method is described 
in detail in Ref. j^j]. The simulation cell consists of a rough bottom in the xy— plane, whose 
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normal is tilted by an angle 9 with respect to the direction of gravitational acceleration g 
in order to induce a flow in the x— direction, as shown in Fig. ^ The bottom is roughened 
by a thin substrate of stationary particles uniformly distributed over the base of the system 
with the same properties as the mobile particles. Periodic boundary conditions are imposed 
in x— and y— directions with a typical base size of 20c? x 20c?, in order to eliminate the 
documented effect of side walls 

The mobile particles flowing above the base are spheres of mass m and diameter c? 
that form Hertzian contacts as well as interacting frictionally. Unless otherwise noted, 
the simulation parameters used for the runs were k n = 2 x 10 5 mg/d (normal contact 
stiffness), /i = 0.5 (Coulomb friction coefficient), k t = (2/7)k n (tangential contact stiff- 
ness), 7 n = 50 \fg~Jd~ (normal viscoelastic constant) and 7t = (tangential viscoelastic 
constant) Q]. In Appendix iBl we also present results with reduced (k n = 2 x 10 3 mg/d) and 
increased (k n = 2 x 10 7 mg/d) stiffness of particles. 

We evolve the system by integrating the appropriate equations of motion with a 
timestep 5t/ro = 10~ 4 , where the characteristic time tq = \J d/g. In order to ensure mea- 
surements in steady state, we equilibrate each system for a long enough time (1 — 2 x 10 7 
St) and check the stationarity of velocity profiles. Our main comparative analysis is applied 
to simulations with a base area of 20c? x 20c? with two different flow heights (h/d ~ 20 with 
9089 particles and 40 with 17889 particles) at different angles 9 in the range between 20.5° 
and 26°. 

The main properties of the runs with typical parameters, as well as the results of the 
correlation analysis for the selected layers in the bulk of these flows, are summarized in 
Table |l] The steady-state volume fraction and velocity profiles of these runs are depicted 
in Fig. El In addition, we simulated and obtained selective results for thin piles, h/d ~ 10, 
with typical (20c? x 20c?) and large (40c? x 40c?) base area (see Appendix [Bj) as well as taller 
piles, h/d ~ 80. 

B. Velocity Correlation Analysis Method 

All correlation analyses are performed at steady state, within thin layers normal to the 
z— axis. The surface correlation analysis was applied to a preset number of particles {N^j^ = 
400 for the typical bases size of 20c! x 20c?, to obtain approximately one layer of particles) with 
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FIG. 2: (a) Volume fraction profiles for all systems listed in Table [I] (b) Velocity profiles for 
h/d ~ 20. (c) Velocity profiles for h/d ~ 40. See also Table [I] for legend information. All profiles 
are obtained with bin size of 0.9d. All data points are connected with lines and the symbols are 
used to mark each data-set. 

the largest z— values. Bulk layers that were analyzed included all particles whose centers 
were located within 0.5<i of the plane located at z = -2i ayer . This yielded A^ ayer ~ 400 — 440 
particles in the layer, depending on 9. The results are not sensitive to variations of N layer or 
bulk layer thickness by 10 — 15%. 

The values for zi ayer , reported in Table HJ are chosen near the center of the pile, in order to 
minimize boundary effects from the top and bottom. In order to ensure that bulk properties 
are observed in these layers, in Fig. El we display the scaling variable I [cf. Eq.®] as a 
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TABLE I: Properties of the simulation runs and associated bulk layers that were analyzed for 
velocity correlations. 

f : From Ref. Q]: The tangential viscoelastic constant = 7n/2. However, h stop is not affected 
by this difference. 

t. : Layer thickness is one diameter, centered around the reported zi ayer . 

function of z— position, along with the position of the layers for which detailed analysis is 
conducted. I is expected to be constant in the bulk Q. 

In the PIV technique used in laboratory experiments, velocities are obtained by mea- 
suring particle displacements between two frames separated by the time At. Similarly, we 
measure the velocity of each particle i G {l..A^ ayer } in a layer from its displacement over a 
characteristic "measurement time" At: 

t£(f) = (r a (t)-r*(f-At))/Af. (5) 

Here a = {x, y, z} is a coordinate label and r l (t) is the position of particle % at time t. 
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FIG. 3: Scaling variable I versus scaled height z/d. 
Due to finite size effects, the instantaneous mean velocity of each layer, 

-^laycr 

Va(t) = ^WMayer, (6) 

i 

fluctuates around the time-averaged velocity profile shown in Fig. I2fb-c). In order to avoid 
introducing spurious correlation effects at large distances, we use the instantaneus mean 
velocity to obtain the velocity fluctuations, rather than the long-time mean velocity: 

<(t)=vi(t)-v a (t). (7) 

Finally, to obtain better statistics, we perform the correlation analysis for « 10 4 con- 
figurations separated by 100 St. Although not explicitly indicated, all velocities v l a have an 
implicit dependence on the measurement time At. 

We compute six two-point velocity correlation functions in the xy— plane for each layer, 
following the method in Ref. Correlations in the (3 = {x, y, z} component of the velocity 
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at a distance r along the a = {x, y} direction are given by 



Here 5(r) = = exp(— ■£^ s ) is a Gaussian function of width w = OAd, as used in Ref. 
and <5 7Q = 1 when 7 = a and otherwise. We do not examine correlations along the 
z— direction since these connect particles with different average velocities and strain rates, 
which could not be statistically averaged in a satisfactorily manner. 

In virtually all cases, a good fit to an exponential decay is obtained. Nonexponential 
behavior was typically related to finite-size effects, poor statistics or dispersion effects, which 
will be discussed in Section ITlI Al A correlation length X a p is associated with each correlation 
function by obtaining the best linear fit to a log-linear plot: 

C at3 (r) = C at3 (0)exp(-r/X a p). (9) 

Note that Pouliquen j<| defined the correlation length as the distance L a/3 for which C a p(r) = 
0.07C ai g(0), which is related to our correlation length through 

L aP « 2.66A a/3 . (10) 

III. RESULTS 

A. Velocity Correlations at the Surface 

Figure |U shows sample vector plots of the velocity field u l (£), projected on the xy— plane, 
with At/ro = 1, for the surface layers of two simulation runs with h/d ~ 20 and incline angle 
9 = 21° and 9 = 23° respectively. Visual inspection of the patterns suggest more correlated 
flow structures for the lower angle run, similar to observations made in Ref. 0. 

For a more quantitative analysis of the velocity correlations, in Fig. El we show all six 
normalized correlation functions (cf. preceding Section) for the pair of runs discussed above. 
From the linear decrease of log [C a p(r)] with distance, we can determine a correlation length 
with each normalized correlation function by fitting to Eq. (|9*jl . 

It is instructive to compare the profiles in Fig. El to those in Fig. El which are obtained 
from exactly the same data set, but by using the time-averaged mean velocity instead of the 



9 



20 



15 



32 

^ 10 



X \\\ 



■H' >- 



\ 1 \ 



20 



15 



32 

^ 10 



10 

x/d 



15 



20 



■ \ ' 



/-/✓ t 



/ ' 



\ 



\ 



I 



\ 



V 



10 

x/d 



15 



20 



FIG. 4: Sample velocity fields in the xy— plane for the surface layers of two runs with h/d — 20 
and At/ro = 1. (a) 6 = 21°, (b) 9 = 23°. 
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FIG. 5: The six normalized correlation functions, C a s(r), for two runs with h/d ~ 20 and 
At /t = 1. (a) 9 = 21° and (b) 9 = 23°. 



instantaneous mean velocity to compute velocity fluctuations. Effectively, this corresponds 
to switching the order of spatial averaging and temporal averaging and is not expected to af- 
fect the result in a large system in steady-state. However, spurious correlations are obtained 
when the latter method is used, which are due to the fluctuations in the instantaneous mean 
velocity shown in Fig. [7| The largest deviations in Fig. Efa) arise for C ax (r), consistent 
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FIG. 6: The six normalized correlation functions, for the same runs shown in FigEl but with 
relative velocities computed by subtracting the time-averaged mean velocity instead of the instan- 
taneous mean velocity. Note the "flattening" of the correlation functions at large distances (see 
text). 
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FIG. 7: Fluctuating components of surface layer velocity for h/d — 20, 6 = 21 c 



with the largest deviations in v x (t) seen in Fig. [7| The resemblance of these to some of the 
correlation functions shown in Ref. j^, which uses the global mean velocity, suggests that 
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the combination of the observed non-exponential tails and the particular definition of the 
correlation length L a p used may have resulted in over-estimation of the reported correlation 
lengths. We typically obtain the least systematic bias in the y— component of the velocity 
and therefore rely mainly on C yy (r) for comparative analysis of runs. 
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FIG. 8: Surface layer correlation length as a function of At for h/d ~ 20 and 9 = 21° (0), 9 = 23° 
(□) and 9 = 25° (O)- 

Although we confirm that the measured correlation is indeed larger for the smaller angle 
run, we find that the velocity correlation lengths depend on the measurement time, At. As 
seen in Fig. |H1 the correlation length initially increases with At, ultimately reaching a max- 
imum value as the mean strain experienced by the particles during the measurement time, 
e = 7 At, exceeds 1. The plateau at large At is understandable: The velocity correlations are 
naturally disappearing as particles start to execute diffusive motion induced by the strain. 
However, the initial increase in X yy makes it difficult to decide the "correct" measurement 
time to use in order to appropriately compare runs with different heights and angles. We 
discuss this issue in more detail in Appendix El and arrive at the conclusion that the results 
should be compared at the fixed values of measurement strain e = j At. 
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B. Depth dependence of velocity correlations 



Since our main interest is in the nature of granular flow in the bulk, we next investigate 
how the velocity correlation length varies as a function of distance from the surface by taking 
advantage of the information available in simulations but usually not in experiments. 
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FIG. 9: (a) Correlation length X yy profiles for angles 23° (squares) and 25° (circles), height 
h/d ~ 20. (b) C yy (r) at the surface (filled circles) and in the bulk (open circles) for 6 = 25° and 
h/d ~ 20. All results are obtained for fixed strain e = 0.1. 



Figure Efa) shows the correlation length profiles as a function of interrogated layer po- 
sition ^ayer for h/d ~ 20 and angles 23° and 25°. The correlation length is slightly smaller 
near the top of the pile compared to the bulk, where it is approximately constant. The 
correlation lengths increase as the angle is decreased towards the angle of repose. This is 
consistent with earlier observations reported on surface correlations |9(. 

Figure HHb) shows typical shapes of correlation functions for h/d ~ 20, 9 = 25°, e = 0.1 
in the bulk compared to the surface. C yy (r) on the surface has more noise than in the bulk, 
which is most probably due to the saltating particles. Fluctuations in the mean velocity 
are dominated by a relatively small number of saltating particles at the surface. Since 
such particles do not exist in the bulk layers, the statistics is better in the bulk layers. 
Furthermore, the rheology of chute flow becomes more robust and better characterized in 
the bulk. Therefore, in order to make a more quantitative connection between the correlation 
length Xyy and deposit function h stO p(0), we will focus on the bulk flow for the remainder of 
this paper. 
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C. Correlations in the bulk 
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FIG. 10: Correlation length \ yv as a function of strain. 



Figure El shows the correlation length X yy as a function of e, for all the runs described 
in Table |U As discussed is Section IIII Al the comparison of correlation lengths at different 
angles should be done at a fixed value of e, for which the rattling motion has largely averaged 
away and decorrelation due to the diffusive motion has not yet set in. For the systems listed 
m Table HI the region 0.1 < e < 0.2 is the most suitable. 

In Fig. ITlTa). X yy , measured at e = 0.1 as a function of /, is displayed. The dependence 
is quite weak, and X yy appears to diverge only logarithmically, if at all, in the I — > 
limit. In contrast, h s t op exhibits stronger dependence and power-law divergence in the same 
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FIG. 11: (a) \ yy versus / and (b) \ yy versus h stop for e = 0.1. 

limit. The interdependence of X yy and h s t op is directly probed in Fig. ITTT b). which reveals 
an approximate relationship \ yy /d ~ ln(h siop /d). This result suggests that the velocity 
correlation length remains much smaller than the pile height for all flowing piles in the 
I — > limit and does not directly influence h stop . In fact, the situation appears to bear an 
interesting similarity to that in molecular glasses, where the diverging viscosity at the glass 
transition is accompanied by a region of cooperative motion that remains quite small and 
appears to show signs of a logarithmic divergence 

To probe smaller values of / one needs to study thick piles at low angles just above the 
angle of repose. However, it has recently been shown by Silbert that these flows exhibit 
flow instabilities. For example, for (h ~ 80) and 9 < 21°, a resonant normal vibration mode 
is present which results in a longitudinal dilation wave in the z— direction. As a result, 
the density in a given layer, as well as the height of the flow, varies periodically in time, 
breaking the time-translation symmetry. 



IV. CONCLUSIONS 



We have presented numerical evidence that the spatial correlations associated with ve- 
locity fluctuations in chute flow decay exponentially with distance and remain short-ranged, 
with the correlation lengths exhibiting at best a logarithmic divergence as the angle of re- 
pose is approached. In order to make meaningful comparisons of correlation lengths, it is 
important to choose measurement times that are inversely proportional to the local strain 
rate. As described in Appendix the analysis of single particle motion reveals two su- 
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perimposed motions: Cage rattling that dominates the instantaneous velocity fluctuations 
but lacks correlation and sterically constrained motion that causes strain-induced diffusion 
and reflects correlations in particle rearragement. A closed-form set of rheological equations 
would need to relate hydrodynamic flow parameters to the microscopic energy dissipation 
in the system, and unraveling this relationship would likely need to take into account the 
quite different dissipation characteristics of these two types of motions. 
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APPENDIX A: DEPENDENCE ON EFFECTIVE STRAIN 

The velocity correlation results depend on measurement time, At, and measurement 
strain e = 7 At. We found it instructive to study the distribution of velocity fluctuations as 
a function of measurement strain e = 7 At. To avoid the effects of boundaries we show here 
the results for bulk layers listed in Table HI Figure El depicts non-dimensionalized mean 
square velocity fluctuations as a function of e: 



where u a is obtained using a measurement time At = e/7. 

The inset in the figure shows that a simple Gaussian distribution is observed in all cases, 
such that the variance completely characterizes the fluctuations. For an individual data set, 
we observe a ballistic-like regime for small e, in which the measured velocities do not depend 
on e, whereas particles exhibit diffusive motion at larger values of e. 

The collapse of all data for large e suggests that the diffusive motion at long times is 
dictated by steric constraints between particles as they pass near each other, such that 
the displacement of each particle depends only on accumulated strain and not the strain 
rate. Hence, the velocity fluctuations arising from this motion do not depend on the scaling 
variable I. 
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FIG. 12: y-component of the mean square velocity fluctuations versus strain, measured in bulk 
layers as characterized in Table [I] Legend is shown in Table [I] Inset shows the probability dis- 
tribution functions for the y-velocity component in the bulk layer for different runs and At's for 
h/d ~ 20. The probability is scaled by maximum value and velocities scaled by variance. 

In the ballistic-like regime the data-sets diverge from each other. In particular, the 
instantaneous velocity fluctuations for small e exhibit anomalous scaling cr^(0) ~ J -1 in the 
limit I — > 0, as shown in Fig. ^] and observed previously [jj. We associate this additional 
motion at short times with the rattling of particles in their temporary cages. A more 
quantitative description of this superimposed grain motion will be discussed elsewhere jl^ . 
This motion is expected to have very weak spatial correlation and average out quickly with 
increasing measurement time, while still dominating the instantaneous velocity fluctuations 
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FIG. 13: Velocity fluctuations versus inertial number / for the systems listed in Table [I] 

and giving rise to their anomalous scaling. Since the overall motion of each particle is the 
superposition of these two motions, it is not surprising that the measured correlation length 
is suppressed at small measurement times. Since we are interested in identifying the nature 
of correlations associated with the sterically hindered motion of the particles, we compare 
correlation lengths for different systems at a fixed value of measurement strain that is large 
enough to have substantially eliminated the effects of cage motion, but small enough that 
the diffusive motion has not caused substantial decorrelation. Note that this effectively 
changes the measurement time that must be used to determine velocities according to the 
local strain rate. 

APPENDIX B: EFFECT OF SYSTEM SIZE AND PARTICLE'S STIFFNESS 

To test the system size effect on the correlation length we equilibrated additional config- 
urations of height h/d ~ 10 at 9 — 23° with two sizes of base area, 20c? x 20c? and 40c? x 40c?. 

Figure IT4T a) shows that the volume faction profiles for the system with b 40c? x 40c? 

is the same as for typical configurations whose bases are 20c? x 20c?. Figure H*4T b) shows the 
velocity correlation functions for the large system compared with the small system size. The 
correlation lengths, obtained from the exponential fit to the first four points, are very similar, 
L w (40 x 40) = 1.295, L yy {2® x 20) = 1.240 , L XX {4Q x 40) = 1.194, L xx {2® x 20) = 1.141, 
(less then 5% difference). The deviation in exponential behavior in the smaller system is 
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FIG. 14: Volume fraction profiles (time-averaged and smoothed) (a) and velocity correlation 
functions (b) for 20d x 20d (□) and 40d x 40d (■) base area with h/d ~ 10, 9 = 23°, At/r = 1. 

larger. The first four points in correlation function are, however, the same. This justifies 
our method of obtaining X yy by fitting the exponential decay function to the first four points 
of correlation data. 
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FIG. 15: (a) Volume fraction profiles for different values of grain stiffness (data smoothing routine 
used in addition to time averaging), h/d ~ 40 and = 23°. (b) Velocity profiles for the same 
configurations shown in (a). 

To study the dependence on spring stiffness k n , we carried out simulations for the tall 
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FIG. 16: (a) Velocity correlation function at e ~ 0.1 and (b) correlation lengths for three values 
of grain stiff nes. h/d ~ 40 and 6 = 23°. 

pile h/d ~ 40 at 9 = 23° with k n = 2 x 10 3 mg/d and k n = 2 x 10 7 mg/d and compared 
the results with earlier k n = 2 x 10 5 mg/d results. Changes of stiffness affect the collision 
time and restitution coefficients and require appropriate adjustments of the timestep and 
viscoelastic parameter !2']- Therefore, for the case of increased stiffness, k n = 2 X 10 7 mg/d, 
the timestep was decreased to 5t/r = 10~ 5 and 7 n was set to 7„ = 500 \fgjd. In the case 
of the soft material, k n = 2 x 10 3 mg/d, the timestep was kept at the same typical value, 
5t/T = 10~ 4 , and the viscoelastic dumping constant was reduced to 7 n = 5 \fg~j~d. 

Figure ITBT a) shows that increasing the stiffness by two orders of magnitude increases 
the value of bulk volume fraction from 0.573 to only 0.575. The stiffer grains flow with 
slightly lower velocities, Fig. ITST b). Decreasing the stiffness by two orders of magnitude has 
a big qualitative effect. The bulk is no longer characterized by a constant volume fraction. 
The grains closer to the bottom are compressed by the weight of the system, leading to a 
decreasing with depth volume fraction profile. Figure ITHT a) shows the effect of the coefficient 
of stiffness on the two-point correlations function C yy {r), measured at At/r = 0.1 for soft 
grains and At/r = 0.25 for typical and stiff grains. This choice of At provides approximately 
equal value of measurement strain, e ~ 0.1. Again, we note a negligibly small difference 
between the results with typical and stiff grains. The results for soft grains indicate a smaller 
correlation length. The difference between correlation length for the three configurations is 
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small and decreases with the increase of e, see Fig. lToT b). In summary, the effect of variations 
in k n is minimal as long as k n > 2 x 10 5 mg/d. 
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